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Abstract. 

In this work we examine the electrostatic screening potential due to a point charge 
Oi ^ located off-centre in a spherical dielectric cavity. This potential is expanded for 

the case in which the dielectric constant e is large, several methods of finding the 
terms in the expansion are investigated, and closed-form expressions are found 
through third order in e along with error bounds. Finally, possible uses of these 
expressions in molecular dynamics simulations of isolated charged molecules is 
discussed. 
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1. Introduction 

Over the last generation, molecular dynamics simulations have become a vital 
theoretical tool in the analysis of the physical interaction of proteins. This has been 
driven in part by the dramatic increase in cheap computer power in the last decade, 
which has grown at almost an exponential rate. 

This growth in computer power has resulted in a similarly dramatic increase in 
the size of the systems studied utilizing this technique. What began as a study of 
modest proteins such as myoglobin has expanded to include complicated systems such 
as structures embedded in cellular membranes and even a tobacco mosaic virus [T]. 

Due to the complexity of these systems, previously ignorable errors due to 
approximations in the model are likely to accumulate, resulting in inaccurate results 
and unstable simulations. Therefore, it is of vital importance to implement as accurate 
a representation as possible, particularly with respect to the long-range interactions 
in the model. Of these, the most problematic is the electrostatic interaction between 
charged elements in the simulation. In addition to generating long-range forces, the 
electrostatic field also polarizes the media in which the simulation is taking place, 
effectively creating more sources for the field in the simulation. It is this aspect of the 
electrostatic interaction that is the most troublesome to implement accurately while 
keeping computational time and expense to a minimum. 

There have been numerous attempts to circumvent the electrostatics problem in 
molecular dynamics models. The classic way of doing this is to place the system in a 
periodic cell and implement Particle Mesh Ewald dynamics [2 to account for the long 
range fields. This model indeed handles the electrostatic problem while keeping the 
system size reasonable. However, it artificially imposes a crystalline structure on the 
system which may not be desirable for some applications. 

If one wishes to investigate an isolated structure, then the options are fairly 
limited. A cutoff on the electrostatic interaction is usually imposed, but this effectively 
isolates portions of the system from one another. These models also suffer from the 
defect that the system in effect becomes finite in size, ignoring a large portion of the 
solvent. Since the solvent — which is usually water — has a large dielectric constant 
(e w 80), it is quite polarizable. Therefore, the field generated by the solvent is very 
sensitive to the background field. Multipole methods ^ historically have had some 
success in dealing with these long range terms. 

We wish to reformulate the approach to the electrostatic interaction in molecular 
dynamics simulations to take into account this sensitive dependence of the system on 
the background field. Our model should have the following properties, 

• It should accurately represent the field of the solvent. 

• It should be relatively inexpensive from a computational point of view. 

• The potential in question should be a solution to Poisson's equation, so that it 
represents a physically possible charge distribution. 

2. Green's Function for the Screening Potential 

If a charge distribution is placed in a dielectric medium that is uniform and infinite 
in extent, the well known result is that the electric potential is reduced by a factor of 
e, the relative permittivity of the dielectric. This reduction is caused by an additional 
"screening potential" due to the polarization induced in the dielectric, which partly 
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cancels the original potential. The problem is more complicated when there are 
dielectric boundaries involved, as in the case of a charge distribution inside a cavity 
within a dielectric. The interaction of charges embedded in a dielectric cavity is a 
surprisingly complicated and rich subject in the study of classical electromagnetic 
theory. Even simple systems fail to have closed-form solutions for the potential. If 
one wishes to construct a realistic model in which charge interacts with a dielectric, 
then some approximation is inevitably necessary. 

To begin building our model, let us assume that a dielectric of uniform relative 
permittivity e fills all of space except for a spherical cavity centred at the origin, 
and that the cavity has unit radius (i.e. all distances are measured in terms of the 
cavity radius). Suppose that there is a charge distribution p(r) within the cavity. The 
potential everywhere in space may be calculated if the Green's function G (r;r') is 
known: 

$(r)= f G{r;r')p{r')d^r'. (1) 

J cavity 

The Green's function satisfies Poisson's equation as a function of r with a unit point 
charge located at r' as the source: 

V2G(r;r') = -47r53(r-r'), (2) 

and must also satisfy the proper boundary conditions on the cavity wall and at 
infinity. This problem is easily solved by the method of images in the limit (e — > oo), 
corresponding to a cavity in a conductor. However, no closed-form expression for 
G (r; r') is known in terms of elementary functions for the case of finite e, except for 
the trivial case in which the charge is located at the centre of the cavity. It is therefore 
necessary to representing the Green's function with an infinite series or other type of 
approximation. 

Let us choose the positive z-axis to pass through the source point r', which we 
are assuming to be within the cavity so that r' < 1. The classic way of tackling this 
type of problem is to write the Green's function as an infinite series in orthogonal 
functions, which for the case at hand will be Legendre polynomials due to the 
azimuthal symmetry of the problem. The potential due to the unit point charge 
alone is 

Gpoint (r; r ) = -j = , , (3) 

|r-r'| y/r^ + r'^ - 2rr' cos 9 

which has a well known expansion in terms of Legendre Polynomials [31|3]. For looking 
at boundary conditions we will be interested in the region near the cavity wall, for 
which r > r' and the expansion is 

oo il 

Gpoint (r; r') = E ;JfT ^' (^os 0) . (4) 

There is also a contribution Gscrcon (r; r') to the Green's function, the screening 
potential, due to the polarized dielectric: 

G (r; r') = Gpoi„t (r; r') + Gscrccn (r; r') . (5) 

Because the effective polarization charge is only on the surface of the dielectric, the 
screening contribution to the Green's function must satisfy Laplace's equation inside 
and outside the cavity, be finite in each of these regions, and satisfy the proper 



Dielectric screening in a spherical cavity 



boundary conditions at the cavity wall r — 1. With these conditions in mind we 
can write 



Gs 



(r;r') 



J2^ir' Pi {cos 0) forr<l, 



;=o 



(6) 



^Air'^^+^^ Pi (cos 9) forr>l. 



l=Q 



The cocfBcients Ai are the same in both sums to ensure continuity of the potential at 
r = 1, which is one of the boundary conditions. 

To find the coefficients Ai, we enforce the other boundary condition, which is 
that the electric displacement be continuous across the boundary. This amounts to a 
condition on the radial derivative of the complete Green's function: 



The result is 



dG 
dr 



Ai = - 



dG 

!_ dr 

{e-l){l + l) a 



(7) 



(8) 



l + e(/ + l) 

Substituting these into © we find our expression for the screening Green's function 
to be 



G(rr',cos6l,e) 



for r < 1, 



Ga 



(r;r')=<^ i .r 



G — , cos 9, e for r > 1, 
r \ r 



where 



G {x, u, e) 



(-i)E 



i + i 



1=0 



l + eil + 1) 



x^Pi {u) 



(9) 



(10) 



for < X < 1 and — 1 < w < 1. This series solution is well known in the literature 
[SllTirsj. We shall refer to the function G [x, u, e) as the "screening function." 

While Equation (jlOp provides a perfectly legitimate expression for the screening 
function, it is fairly limited to its usefulness in numerical simulations of charges within 
the cavity. The rate of convergence depends on where the charge is within the cavity, 
and is rather slow unless the charge is close to the origin. However, we can look for 
methods of summing the series with hopes of obtaining an expression that may be 
more useful in calculations. The coefficient within the sum in (1101) can be written as 



/ + 1 



1 



+ 



1 



1 



; + e(; + l) e+l ' (e + lf l + e/{e+l)' 
separating the sum into two parts: 



G (a;, u, e) = — 



1 



£ + 1 



oo ^ oo 



1 



l + e/{e + l) 



x'Pi {u) 



.1=0 1=0 

The first is just the generating function for the Legendre polynomials [9], 

1 



J2x'Pi{u) 



1=0 



y/l — 2ux + x"^ 



(11) 



(12) 



(13) 
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The other sum can be evahiated in several ways, one of the simplest of which is to 
start with the generating function again, 

multiply both sides by t'^/('^+^)^^^ and integrate: 

Vf^ -2tu + i Jo r-'n 



1= 



,^„« + e/(e+i: 



E 



Piiu). (15) 



This allows us to express the Green's function in terms of an integral|: 



G{x,u,e) = — 

e + 1 



1 a;-^/(^+i) r t-i/(^+i)di 



, , , , . (16) 

Vl - 2ux + x^ 1 + e Jo Vl-2ut + t^\ 

Expression p^ lends itself to the following physical interpretation: First of all, in 

the limit e — > oo, where the dielectric becomes a conductor, only the first term in the 

square brackets survives. For the solution inside the cavity, this corresponds to the 

familiar image point charge outside the cavity, and for the solution outside the cavity, 

this corresponds to putting an image charge at the location of the original charge, 

neutralizing the original charge and "grounding" the conductor. Now for finite e, the 

image charge is reduced in magnitude by the factor (e — l)/(e + 1) so that there is a 

nonzero potential outside the cavity, but this point charge is not enough to satisfy the 

boundary conditions; an additional line image charge is needed, stretching from the 

image point charge to infinity (for the inner solution) or to the origin (for the outer 

solution) . 

Figure [1] shows a plot of the screening function as a function of the variables 

Xi = xcos9 and X2 = xsinO in the circular region defined by < a: < 1- For points 

inside the cavity, x — rr' as indicated in ([9]), so xi and X2 are proportional to the 

spatial coordinates parallel and perpendicular to the z-axis, respectively. For a given 

location r' of the source charge, the cavity boundary corresponds to a; = r' which is at 

most one, so a plot of the potential inside the cavity may be visualized by truncating 

the plot in figure [1] to a smaller circular region of radius r' , and then expanding the 

plot to fill a region of radius one. The screening function becomes infinite at the point 

(a; = 1,M = 1), corresponding to the image point charge referred to in the previous 

paragraph, but this singularity only shows up in the physical solution when r' = 1, 

i.e. when the source charge is at the cavity wall. 

3. Expansion for Large Values of the Dielectric Constant 

Equation p6p provides a concise, exact expression for the screening function, which 
can be evaluated numerically and also in terms of Appell hypergeometric functions 
[lOj . However, it is of limited use in actual simulations of molecules in which great 
numbers of these expressions would need to be evaluated. It is therefore useful to look 
for approximations allowing the use of simpler functions. One such approximation is 
to recognise that, as we have mentioned, the dielectric constant of water is quite large, 
so that an expansion good for large values of e would be useful. 

f This equation is equivalent to the one appearing in references [3 18] . 



Dielectric screening in a spherical cavity 



G(x,u) 




Figure 1. A plot of the screening function for e = 80, produced by numerically 
evaluating and plotting JTBJ using the standard packages in Mathematica. The 
axis coordinates in this and subsequent surface plots were chosen to be xi = xu = 
X cos 9 and X2 = ix\/l — u'^ = ±x sin 9. 



Based on what we have aheady derived here, there are two ways in which we can 
obtain a series expansion good for large e. One way is to go back to definition p^ of 
Gn (a;, u) and write it in the form 

-1 

1 - T^, X^Pl (U) 



£-1 



e- 1 



1=0 

oo oo 



ie+l){l-l)_ 



e + 1 ^^^^ ie+lYHl + l) 



Reversing the order of summation gives us a series in powers of l/(e + 1)": 

G {x,u, e) 

where 



e-1 / Gi{x,u) G2ix,u) Gsix^u) 

Go[x,u) + 



e + 1 



e+1 



(e+lf {e+lf 



1=0 ^ ' 



(17) 



(18) 



(19) 



The other way to get a series is to expand the integral expression P^ . noting that 



n-l 



1 1) 1 ■^-^ in 

1 ^ x^(k-\ 



ln""^(x/t) 



X e+1 a;-^^ (fc-l)!(e+l) 



— ■ (20) 



This gives the series in p^ again, where this time the coefficients are given by 

1 



and 



Go (x,u) = 



G„ (x,u) 



\J\ — lux + x^ 

1 r ln"-i(x/i) 



1 



{n-\)\ xJo VI - 2ut + t^ 



:dt 



for n> 1. 



(21) 



(22) 
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Expression (pij) is clearly equivalent to series expression P^ for n = 0, since the 
former is the generating function for the Legendre polynomials. It is also possible to 
show directly that series expression ([T^ and integral expression ([^^ are equivalent 
for n > 0. Our purpose in deriving these two expressions for Gn {x, u) is that they 
complement each other: The series is useful for deriving general properties of these 
coefficients, while the integral is useful for doing actual calculations of them. 

To get an idea of what sort of functions may be involved in finding the coefficients 
Gn {x, u), let us use the series expression to evaluate them in the special case u = ±1, 
corresponding to points on the z-axis. Since Pi (±1) — (±1)', equation P^ gives 

G„(x,±l) = 5](±l)'-^ = ±-Li„[±a:], (23) 

where Li^ is the polylogarithm function [TTJ [T^l [13] : 

oo ^ 

LinW-Elr- (24) 

fc=i 

As can be easily seen from this definition, the polylogarithms satisfy a recursion 
relation, 

Li„N= r^^ili-iMdi. (25) 

Jo t 

The polylogarithm for n = is easily found from definition (j24p : 

LioH = -^, (26) 

1 ~ X 

and those for n — 1, 2, 3, . . . may be obtained from this function by successive 

integrations. As a consequence, 

Lii[2;] = -ln(l-x), (27) 

while the dilogarithm Li2 [x] , the trilogarithm Lis [^] i ^iid higher-order polylogarithms 
cannot be expressed in terms of elementary functions. For any positive integer n, 
Li„ [x] has a real value for — oo < x < 1, but has a branch cut in the complex plane 
running along the positive real axis from x = 1 to oo, across which the function has a 
discontinuous imaginary part. 

Like the polylogarithms, the coefficients G„ {x, u) also have a similar recursion 
relation which follows immediately from series expression (|19p : 

Gnix,u)^- I Gn-i{t,u)dt, (28) 

and in fact the recursion relation for xGn {x,u) is the same as for Li„ [x]. Equation 
(|28p serves as an alternate way to calculate G„ (x,u) by starting with Gq {x,u) as 
given by ([2T|) . 

To calculate G„ {x, u) for specific values of n, we can either use integral expression 
(|22p or recursion relation (P5|) : both methods seem to lead to about the same degree 
of complexity. In both methods we have found it useful to make the change of variable 



l + t-y/l-2ut + t^ w ( \ + u 

w = -— , or t^- 1 ^— w . 29 

1 + u 1 — w \ 2 

In addition, we define the following variables for the sake of convenience. 



1 + X — y/l — 2ux + X^ 1 — M 

1 + M 1 + M 



Dielectric screening in a spherical cavity 



For physical values of x and u, the quantities p and q lie in the range < p < 1 
and < (7 < cx)l§| The first-order coefRcient Gi {x,u) is then easily evaluated using 
either the integral expression or the recursion relation; both methods lead to the same 
integral: 

^ , , 1 /■" dt 1 fP dw 

Gi [x,u) 



X Jo 
1 



Hi-p) 



2ut + i2 
1 



1 



(31) 



■ Lii [p] , 

X X 

The higher-order coefficients involve higher-order polylogarithms and become 
increasingly complicated. A pattern that emerges is that for n > 1, xGn {x, u) can be 
expressed entirely in terms of logarithms and polylogarithms of rational functions of 
p and q. Using transformation ([29ll on integral ([22|) . along with ([30|) to eliminate x in 
the integrand, we obtain 



xGn ix,u) = 



1 



{n-iy. 



w 



bil-w) 



w{l 



w) 



dw 



1 



(32) 



where b = p(l 
relation, 



q — p)/{l — p). Also, from ((28)) we obtain a transformed recursion 



xGn {x,u) ^ - 

X 



P /I 

tGn-l{t,u) - + 



1 



1 



1 — w 



dw, (33) 
/O \uj ±—w 1 + q — wJ 

where tGn-i (t, u) is assumed to be written in terms of w and q. Armed with either of 
these equations we can evaluate the second-order coefficient without much difficulty. 
For example, using ([55)) we have 



XG2 {x, u) 



Pln(l 



w) 



dw 



P ln(l - w) 



P ln(l - w) 



dw. 



/o uj Jo J^ — "7 jQ 1 + q—w 

The first integral is just Li2 [p] and the second one is easily evaluated as -^ In (1 
The third integral can be found by making the change of variable w' — q/(l + q 
the result is — Li2 [w'] — ln(w' /q) evaluated at the endpoints. The final result is 



(34) 

-P)- 
-w): 



G2 {x,u) = - <^ Li2 [p] 

X 



Lio 



i + q 



Lio 



1 



+ ^ln^{l-p) + lln^{l + q) 



p + q 
1 



ln2(l 



-p + q) 



(35) 



The third-order coefficient is considerably more complicated, whether we use ()32|) or 
(|55)). We evaluated it using ([5^ . enfisting the aid of Mathematica [I4 to find and 
manipulate the large number of terms and help simplify the expression. In the process 
we also made use of a number of dilogarithm and trilogarithm identities [12[ I13| . Our 
result is 

G^{x,u) = -{2hh[p]-U. 

X 



Li.s 



-Li3 


P 

1 + 


q. 


+ L 


3 [1 - p] - Li3 


'1-p + q 

[ 1 + 9 \ 




[ ^"^ 1 

l-p + q_ 


+ Li3 


1 


-Li3 


q 

l + q\ 


+ Li3 


U + q\ 



1 -p + q 



§ The ranges of these variables are interconnected by the fact that, for a given vahie of q, the 
maximum value of p (corresponding to re = 1) is 1 + g — \/q{l + q), which is 1 in the limit m -+ ±1 
but less than 1 otherwise. 
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Lis 



pq 



6 



■In^ 



l-p 
1- 



In 



p + q 



jl-p + q? 
il-p){l + qy 



(1- 
lnpln(l 



P)(l 



^9) 
ln(l 



2 1 + 



q)lii 



■In^ 



Li2 



1 



1 + 9 



(36) 



9 



(1 



1 



In 



-P)(l- 
l-p 



q 



il~p)il + q) 



One of the complications in writing down the coefficients past first order is that 
there are many different-looking ways to express each G„ {x,u), due to the numerous 
identities satisfied by the poly logarithms. In an attempt to be somewhat systematic in 
our expressions, and also to facilitate the study of their behaviour as well as repeated 
calculations that would occur in plotting or numerical simulations, we have used these 
identities where necessary to express the results in a form where the arguments of all 
polylogarithms are between and 1 inclusive. 

As a check on these results, we can look at their behaviour in the limit u ^ 1, 
which corresponds to 

p ^ X, g ^ 0. (37) 

In ([55)1 , the first term in the curly brackets becomes Li2 [x] while the other terms either 
become zero or cancel in pairs, so that G2 (x, 1) = Li2 [x] jx as expected from (p3)) . In 
P6p . the first two terms give Lis \x\ while the rest becomes zero, giving the expected 
result G3 (x, 1) — Li3 \x\ jx. The behaviour in the limit u ^ — 1, while correct, is less 
straightforward; it corresponds to 



V 



and in this limit wc obtain 



00, 



G2(X,-1) = - Li2 



1 



+ -ln2(l + a;) 



(38) 



(39) 



which reduces to the expected — Li2 [— a;] jx by way of one of the dilog identities. The 
behaviour of the third-order coefficient is similar. 

It should be pointed out that despite the factor of 1/x occurring in these 
expressions, G„ (x, u) is not singular at a; = 0, because the rest of the expression 
becomes zero there. In fact, 

G„(0,u) = l, (40) 

as can be seen from series expression ([T9| . 

Figure [H shows three-dimensional plots of Go(a;,cos0) through Gzix^costH) as 
fimctions of x\ = xcosQ and X2 = a; sin 6*, similarly to figure [1] The zeroth- 
order coefficient has an inverse first-power singularity at the point (x = 1, u = 1) 
corresponding to the image point charge as per the discussion following equation 
(|16p. and the first-order coefficient has a logarithmic infinity at that point, while the 
rest of the coefficients are finite everywhere in the region. Despite the increasing 
complexity of the expressions as the order becomes higher, their actual behaviour 
becomes increasingly simple; the plot of the third-order term is comparatively flat. We 
can in fact obtain bounds on G„ (x, u) in general by noting from integral expression 
(|22|1 that for any given x, the integral will have a maximum value when u — \ and 
a minimum value when m = — 1. It therefore follows from (I23p that for given x, 
the maximum value is Li„ [x] j x and the minimum value is — Li„ [— x] j x. Since for 
< X < 1, Li„ [x] /x and — Li„ [— x] jx are strictly increasing and decreasing functions 
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Figure 2. The coefficient functions Go through G3. 

of X respectivelylill it follows that (for n > 1) the maximum value over all x and u is 
Li„ [1] — C("), and the minimum value is — Li„ [—1] = (1 — 2^"+-*^) C(")j where ((n) is 
the Riemann zeta function: 

C{n) = ^ + ^ + ^ + --- . (41) 



V 



3" 



It follows that 



(l-2-"+i)C(n)<G„(2:,u)<C(r^) 



(42) 



for all X and u in the physical region, and for n > 1. The n = case is easily treated; 
the lower bound is Go (I7 — 1) = \- The results are numerically, 

0.500000 < Gq(x,u) < 00 

0.693147 <Gi(x,m) < cx) 

0.822467 < G2 (x, u) < 1.644934 (43) 

0.901543 < G3 {x,u) < 1.202057. 

As n becomes large, the lower and upper bounds both approach 1, leading to an 
increasingly flat plot. The coefficient Gn{x,u) could then be approximated by a 
simple polynomial, or even a constant, greatly reducing time in computationally 
intensive problems. However, for problems occurring in practice with large values 
of the dielectric constant, it may not be necessary to keep very many terms anyway. 
We will look at the errors due to truncating the series in section [5l 

4. Electric Fields 



While the potential is useful for calculating the potential energy of the system, one 
also needs to be able to calculate the force of the various objects upon one another 
so that the system may evolve from one time step to the next. This in turn requires 
knowing the electric field due to the charges and the dielectric. Corresponding to the 

II This can be seen by looking at their derivatives, using the power series in I I23I I. 
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Green's function G (r; r'), which is the potential at r due to a unit point charge at r' 
in the presence of the dielectric, let us define E (r; r') be the electric field at r due to a 
unit point charge at r' in the presence of the dielectric. The electric field is determined 
from the Green's function by 

E(r;r') = -VrG(r;r'), (44) 

where Vr means the gradient with respect to the vector r, keeping r' constant. Like 
the potential, the electric field separates into the field to the point charge alone plus 
the "screening field" due to the polarized dielectric: 

E (r; r') = Epomt (r; r') + E.erccn (r; r') . (45) 

The field Epoint (r; r') is the familiar Coulomb field, 

F (r- r-'\ - ^^^' - (r-r'cos6>)e,, + (/sin6>)eg 

I^pointir,rj |j._j.,|3- (^2+^/2_2rr'cOS0)3/2 ' ^ ' 

where e^ and eg are unit vectors in the direction of increasing r and 9. In ^ we 
expressed the screening part of the Green's function inside and outside of the cavity 
in terms of a single function G{x,u,e), which can be thought of as a function of 
a three-dimensional coordinate x with spherical coordinates {x,9,(f)). Let us define 
E (x, u, e) to be the negative of the gradient of this function: 

d 1 d 



E {x, M, e) = - VxG (x, u, e) = - ( e^, 




= l-e - 

\ ^ dx " X 

Then from ([9]), the screening field can be written as 

{r'E(rV, COS0, e) for r < 1, 

r' -/r' \ I -fr' \ (48) 

— E I — ,cos6', e I H — j G I —, cos 0,6 I e^ for r > 1. 

The expansion of G(a;,u, e) in powers of 1/(1 + e) leads to a corresponding 
expansion of E (x, u. e): 



^{x,u,e) = --^\ Eo (x, u) + \\ \ ' + /^ [J + /^ .; + . . . , (49) 



e - 1 / Ei(x,m) E2(x,m) E3(x,u ) 

Eo (x, u) -\ ; 5" H 

£ + 1 V ^+1 (e + 1) (e+1) 

where 

E„(x,u) = -VxG„(x,m), (50) 

The zero-order term is found in a straightforward manner, and is just the Coulomb 
electric field due to a unit point charge located on the x-axis at unit distance from 
the origin: 



[x - u) e^ + ^/l^^ Bg 
^°(^'")= (l-2ux-t-x^)3/^ ■ ^'^) 

For the higher orders, we can avoid having to take some of the derivatives of G„ (x, u) 
explicitly by noticing that if we differentiate both sides of recursion relation (|28p with 
respect to x, we obtain 

BC 1 

^--(G„-i-G„) forn>l. (52) 
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The derivative with respect to u is not so easily found, but even it can be siniphfied. 
For n > 1, we found that our expressions for G„ {x, u) were of the form 1/x times a 
function consisting of logarithms and polylogarithms of p and q. Therefore, for n > 1 
it is convenient to write 



dGn 



1 d{xGn) 

X 



dp djxGn) _|_ dq_ djxGn) 
du 



dp 



du dq 



du X du 
On the other hand we also know that, again from the recursion relation, 

d{xGn) dp djxGn) _^ dq djxGn) 
dx 



Gn~l — 



(53) 



(54) 



dx dx dp dx dq 

From definitions (j30p of p and q, the partial derivatives of them with respect to x and 
p are 



dp 


(1- 


pfi^ + Q) 


dp 
du 


P\l' 


p){^ + q) 


dx 
dq 
dx 


(1- 
= 0, 


~pY + q ' 


2[(1^ 
dq _ 
du 


- p? + q] 

[l + q? 
2 



(55) 



Since dq/dx — 0, equation ([M]) allows us to express d{xGn)/dx in terms of Gn^i, 
which can then be used in (|53p . giving 

dq d{xG„ 
du 



dG„ 



dt 



dp/du 



dp/dx 
p' 



2(1 -p) 



Gn-l 
Gn-1 



dq 

[l + qf d{xGn) 



(56) 



Using these results in (1171) along with the fact that vT" 
for the n^'^-order coefficient in the electric field 



,2 - 



= 2y/g/(l + q), we obtain 



E„(a;,u) = -(G, 

X 



G„ 



\ , Vq 

_i)e, + — 



p 



Gn-1-il + q) 



d{xGn 



(57) 



(l-p)(l + g) 

so that we only need to take derivatives with respect to q. It should be noted that, like 
Gn {x,u), E„ {x,u) is not singular at a; = despite the inverse powers of x, as they 
multiply factors which become zero there. In simplifying expressions for the potential 
and electric field it is sometimes convenient to use the expressions for x and Gq in 
terms of p and q: 

_ p{^-p + q) 



il-p){l + qy 



Go {x,u) = 



1 



VT^^i 



(i-p)(i + g) 

(l-p)2 + g ■ 



(58) 



(59) 



As an example in finding a higher-order coefficient in the electric field, let us 
evaluate Ei{x,u). From (pij) . xGi = — ln(l — p) which is independent of q, so no 
derivatives have to be taken at all in evaluating ((57)) . After some simplification we 
obtain 

p^Vq 



El {x,u) = 2 



p{l-p + q) 

(l-p)2+(? 



+ \n{l-p) 



1 



.x^ (1 — p)2 + q 



(60) 



The (1 — p)2 + <z in the denominator will be recognized as the inverse first-power 
singularity characteristic of Go {x,u). As with the potential, the electric field is finite 
at a; = despite the inverse powers of x. 
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5. Analysis of Errors 

In implementing our series expansion (|18p for any particular problem, we need to know 
when we can when we can truncate the series for a given desired degree of accuracy 
(particularly considering how increasingly complicated the series coefficients become 
with each subsequent order). Suppose that we have evaluated the series from n = 
to n = N. The error e^ due to ignoring the terms from A^ + 1 to infinity is 



-N 



e- 1 v^ Gn {x, u) 

e+1 ^ (e + l)» ■ ^ ' 

n=N+l ^ ' 

The terms in this series are all positive, and the maximum value of G„ {x, u) is C("-) 
from (gll), so 

n=N+l ^ ' 

This sum can be evaluated numerically for any given e and N. For example, using 
e = 80 we find 

ei < 2.468 X 10"^ 

62 < 2.231 X 10"*^ (63) 

eg < 2.482 X 10"^ 

(64) 

We can also find an upper bound for the fractional error e^/|G (x, u, e) | by noting 
that the true value of \G {x, u, e) \ is greater than the iV*''-order calculated value, and 
therefore the fractional error will be less than that found by dividing the error by the 
calculated value. Also, a lower bound over all x and u for the calculated value can 
be found using the lower bounds for G„ (a;,^) in (|15|) . For example, suppose that we 
have kept only the n — and n — 1 terms, again assuming e — 80. The maximum 
error is e^ in (|63p . while the minimum value of \G {x, u, e) | is at least 

79 / 0.6931\ _ 



\G {x, u, e) Uin = — (0.5000 + -^^1 = 0.4961 (65) 

and so an upper bound for the fractional error for any x and u is 2.468 x 10^*/0.4961 = 
4.974 X 10~^. Similarly, for a calculation through second order the fractional error 
will be less than 4.496 x 10~^, and for a calculation through third order, less than 
5.001 X 10~*. Of course, the fractional errors will be much smaller than these bounds 
in regions where G is large, even though the error itself do not vary much over the 
whole region. 

6. Use in Molecular Dynamics Simulations 

Having obtained our expression for the potential to third order in (1 + e) , we now 
turn to the question of feasibility. That is, to what degree is this calculation useful to 
those doing large-scale simulations of biological systems? 

A cursory objection may be that the computation of special functions as required 
by this potential for each of the charged objects in a simulation would be intractable, 
and that it would lead to an inexorable increase in computation time. While certainly 
true, the problem is not as severe as it first seems. 
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First, the computation of polylogarithnis is very well understood. The Taylor 
series converges rapidly for arguments in the range < x < 1. Moreover, a partial 
fraction technique exists J15l that converges extremely rapidly. For instance, an 
evaluation of Li2 [^] converges to within 1 x 10~^ of the actual value utilizing only 10 
divisions. 

Second, if the system is of a suitable size, some fields can be calculated using 
the total charges of extended structures, such as the residues comprising the proteins, 
rather than the constituent atoms. Since each term in the potential is separately a 
solution to Laplace's equation, they each comprise a complete solution. For instance, 
one could approach a large system by calculating the atomic charges exactly to 
O [(1 + e)^^], while using the G2{x,u) term to calculate screening by using total 
charge of the protein residues and treating them as point charges. 

7. Summary 

For the potential of a point charge in a spherical dielectric cavity, the authors have 
computed integral expressions that lead to a calculable expansion in inverse powers 
of the relative permittivity. The main feature of this expansion is that the truncated 
series yields a potential that is accurate over the physical region while remaining a 
solution to Poisson's equation at each order. For the case of water, truncating the series 
to second order leads to a fractional error of no more than 4.5 x 10~^, and truncating 
the series to third order leads to a fractional error of no more than 5.0 x 10~^. 

While much of the theoretical groundwork in this article is complete, the 
feasibility of this model in an actual simulation environment has not been studied. 
The polylogarithms introduced to describe the series coefficients G„ {x, u) may be 
computed rapidly and precisely using the partial fraction technique, but it remains 
to be seen whether the accuracy gains surmount the time lag introduced by having 
additional terms in the potential. A systematic study of this would need to be done 
by altering one of the existing molecular dynamics packages such as CHARMM [16] 
or NAMD [Hj. 

Appendix A. Contour Representations 

In addition to the methods that we have presented for finding the screening function 
and the coefficients in the expansion for large e, there are other methods using complex- 
variable techniques. For example, the second sum in (fT^ can be found by looking at 
the special case u — \, which can be expressed as a hypergeometric function 9 : 



'-^^^x^^Kf. 



'+^t,l 



.1;1 



r ' e + r 



(A.l) 



e+l 

The Legendre polynomials can then be brought in by using a contour-integral 
expression for them, 

Pi{u) = —i — , A.2 

2tt\ J y/l- 2uz + z^ z 

where the integration is counterclockwise along any path that encircles the origin but 
does not encircle the branch points of the square root (e.g. a circle of radius less than 
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1). The sum with the Legendre polynomials can then be expressed in terms of a 
contour integral: 



E 



0^ + 



-Pi (") 



e+1 



1 

27ri 



E' 



1 



{z/xY 



dz 



1=0 i + 
1 + e 1 



£ + 1 

2F1 



Vl - 2uz + z2 z 



^'i;i 



£ + 1 ' Z 



e 27ri / Vl - 2uz + z"^ 

The screening function then becomes 

11 /• 2-Fi 



dz 



G (a;,u, e) 



1 



1 



1 



\/l — 2ua; + .T^ 



^'i;i 



Vl - 2uz + z2 



dz 

z 



(A.3) 



(A.4) 



A similar technique can be used in expressing the n'^-order coefficient G„ (x, w), 
using series expression (j23p for the case m = 1 in combination with contour integral 
(US, giving 

"a;" 



G„(x,u) = -— ^ 



Lir, 



: dz. 



(A.5) 



X 27ri / Vl - 2zM + z2 

All of the expressions for the G„ listed in the main body of the text may be obtained 
by directly evaluating this contour expression. In particular, it is easy to derive the 
recursion relation ([28|l by direct differentiation of (jA.5|l . 
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